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We study the lattice 0(2N) Gross-Neveu model with Wilson fermions in the fermion loop for- 
mulation. Employing a worm algorithm for an open fermionic string, we simulate fluctuating 
topological boundary conditions and use them to tune the system to the critical point. We show 
how the worm algorithm can be extended to sample correlation functions of bound states involv- 
ing an arbitrary number of Majorana fermions and present first results. 
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1. Introduction 

It is well known that fermionic degrees of freedom are difficult to simulate on the lattice due 
to their Grassmannian nature. The fermion matrix determinant, obtained after integrating out the 
fermion fields, remains the major bottleneck of Monte-Carlo simulations, since it is highly non- 
local and connects all the degrees of freedom. Each Monte-Carlo update step changes the fermion 
matrix and requires a new calculation of the determinant. On the other hand, if an algorithm 
can update the determinant by a local procedure, then a significant gain in efficiency might be 
achievable. However, simulations of fermions usually suffer from critical slowing down when the 
fermion correlation length diverges, i.e. when the fermions become massless, and it is by no means 
obvious whether a local algorithm can be constructed to circumvent this problem. 

An interesting attempt to deal with the problems related to simulating fermions on the lattice 
has recently been suggested in [1, 2] in the context of the 0(2N) Gross-Neveu (GN) model. It is 
based on reformulating the model in terms of closed fermion loops [3, 4, 5]. Introducing an open 
fermionic string, an algorithm can be devised for which critical slowing is essentially absent. It 
follows the spirit of Prokof 'ev and Svistunov's worm algorithm [6] and makes use of the fact that 
a global update of the closed fermion loops can be obtained by locally updating the open fermionic 
string. The open string corresponds to the insertion of a Majorana fermion pair and directly samples 
the two-point correlation function. In this way the configurations are updated on all length scales 
up to the correlation length, and this eventually guarantees the absence of critical slowing down. 
Moreover, the algorithm also allows simulations directly in the massless limit [1,2] and provides 
direct access to the critical point via ratios of partition functions [7] . 

The algorithm has been successfully applied to simulate free fermions in two dimensions, but 
it has proven equally successfull also in its application to strongly interacting fermions, e.g. in 
the Schwinger model in the strong coupling limit [1,2] and in supersymmetric quantum mechanics 
[8]. Here we report on the application of the worm algorithm to another two-dimensional system of 
interacting fermions - the 0(2N) Gross-Neveu model. We demonstrate how a simple modification 
of the worm algorithm can be used to measure the fermion bound state spectrum and we present 
preliminary results for N = 1 , in which case the GN model corresponds to the Thirring model. 

2. Fermion loop formulation of the Gross Neveu model 

We consider the two-dimensional 0(2N) -symmetric GN model [9] described by the Lagrangian 

& = E Wiiy^ + m) ¥i— y(E wmf • (2- d 

i= 1 i= 1 

This is a relativistic quantum field theory of N self-interacting Dirac fermion fields. For N = I the 
0(2) GN model is equivalent to the massive Thirring model as (yy) 2 = \(WYhY)(WYhY)- It has 
a pseudoscalar bosonic fermion- antifermion bound state [10] and is especially interesting due to its 
equivalence to the sine-Gordon model [11], in which the boson is the fundamental particle and the 
fermion emerges as a soliton solution. 

To make the model amenable for the fermion loop algorithm in [1], we decompose the com- 
plex Dirac fermion fields \jfj and Xjfj into real Majorana fields £2/ and %2j+\ according to = 
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(i§2; ■ + /<§2y+i)/v / 2 and Xj/j = (^ 2 j~^2j+l)^- The Lagrangian density can then be written as 

1 2N - „2 / 2 _ \ 2 

Note that the 0(2N) symmetry is now obvious, as (2.2) is invariant under rotations of the Majorana 
fields in flavour space. Discretising the action with Wilson fermions gives 

1 2N _ p 2 / 2N _ \ 2 2N _ 

s= «EE^*W6W-|-E E«i(x)6(x) -EEfotoWSKx+A)), (2-3) 

Z * i=l 8 * \i=l / i=ljc,jU 

where q> = (2 + m) andP(± i u) = i(l^y^). When expanding the Boltzmann factor in the partition 
function to all orders, terms quadratic or higher order in each field component vanish due to the 
nilpotency of the Grassmann fields. Restricting now to N = I for simplicity and keeping only 
non-trivial terms, the partition function can be written as 

z = / ^nf 1 -! (£?/(*)&(*)) +^^«iW«iW?2W&w) 

* V v-i / / (2 4) 

x,n i=l K ' 

The integration measure is saturated site by site by combinations and it is straightforward to 
identify the non- vanishing contributions to the partition function. They are characterised by the fact 
that for a particular Majorana flavour either two adjacent hopping terms and no monomer terms, 
or one of the monomer terms, but no hopping terms are present at a given site. This results in the 
constraint that for each flavour only closed, non-intersecting fermion loops survive the Grassmann 
integration. Furthermore, the loops are non-backtracking due to the orthogonality of the projectors, 
viz. P(+jj.)P(— ju) = 0. Thus, the partition function is a sum over all possible combinations of two 
different species of loops (corresponding to the two Majorana flavours). The generalisation to an 
arbitrary number of Af Dirac fields is straightforward - the number of species of loops involved is 
simply equal to 2N, the number of Majorana fermions. 

3. Worm algorithm for Majorana fermions 

In order to generate configurations of closed loops we employ a variant of the algorithm of 
Prokof 'ev and Svistunov [6]. Here we explain the main ideas of the open fermionic string ("worm") 
algorithm in a few schematic steps and point out the modifications we have introduced to increase 
efficiency. Further details can be found in [6, 1]. 

A peculiar feature of the worm algorithm is that the fermion correlation function is measured 
during the update procedure. This is due to the fact that the insertion of the open fermionic string, 
which is used to update the loop configuration, corresponds to the insertion of a pair of Majorana 
fermions £i(x)%i(y) of flavour i at positions x and y, respectively. In the path integral formalism 
this is equivalent to the correlation function 

G,(x,y) = Ux%{y)) = y®$ $i{x%{y)e- S ■ (3.1) 
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Figure 1: Update moves for the break-up/reconnect step. Note that the configuration in the middle has no 
physical interpretation and hence gives no contribution to the correlation function. 

The algorithm now proceeds by locally updating the ends of the open string using a simple Metropo- 
lis procedure according to the weights of the corresponding two-point function. The main steps are 
as follows: 

• Relocation step: Given a closed loop configuration, choose a fermion flavour i and a lattice 
site x at random and place the head 5 and the tail 5 of the worm on this site with a probability 
given by the weight ratio of the loop configurations before and after the step. If accepted, 
it gives a contribution to Gi(x,x) unless a loop of species i passes through x in the original 
closed loop configuration, in which case the new configuration has no physical interpretation 
and gives no contribution to G,(x,x). 

• Move step: Choose a direction pt at random, and move the head of the worm to site y = x + jX. 
Add or delete a fermion bond between x and x + {X depending on whether the bond is empty 
or occupied. The resulting configuration gives a contribution to d{x,x + pt). 

• Break-up/reconnect step: In case the new site x + /t is already occupied by a fermion loop of 
flavour /, we still allow the move, although the corresponding configuration (cf. middle plot in 
figure 1 ) is forbidden by the Pauli exclusion principle. Consequently, it does not contribute to 
the correlation function, but induces transitions between allowed configurations as indicated 
in the figure. 

• Removal step: Once the head of the worm reaches its starting position, i.e. head and tail meet 
again at site x, we may propose to remove head and tail. If accepted, we have a closed loop 
configuration contributing to the partition function. 

We emphasise that the break-up/reconnect step is crucial for the algorithm to work efficiently, 
since it allows the loops to be opened and restructured. Especially close to the critical point, where 
loops proliferate, this step gives the worm much more freedom to update the configurations. 

4. Topological and fermionic boundary conditions 

On a finite lattice with a periodic torus geometry, fermion loops can wind around the lattice 
and the loop configurations can hence be categorised into different homotopy classes depending 
on the number of loop windings. For each Majorana fermion <^,- a two dimensional vector /; = 
(l x ,h)i is assigned to each configuration to account for the windings in space and time direction, 
respectively. The components of /; are either or 1 corresponding to an overall even or odd number 
of loop windings, respectively, in the corresponding direction. Configurations with different /, 
contribute to separate partition functions Zr ? with fixed topological boundary conditions (b.c). 
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Figure 2: The figure on the left shows the ratio Z 0100 /Z for the Thirring model at g = 0.9 as a function of 
the bare mass am on square lattices of size L = 8, 16,32,64, 128 and 256. The inset shows the finite size 
scaling. The right plot summarises our results for the critical mass am c as a function of the coupling g. 



Since the open fermion string of the worm algorithm tunnels between the configurations in the 
various homotopy classes, it samples all the partition functions Zj j . More importantly, the 
configurations in all homotopy classes are sampled with positive weights relative to each other, i.e. 
the partition function Z = Y,f j 2 Zjr f corresponds to one with fluctuating topological b.c, but 
unspecified fermionic b.c. As a consequence, the b.c. for the fermions can be chosen at the end of 
the simulation and all possible fermionic b.c. can be studied a posteriori. 

In order to encode the fermionic b.c. we introduce a two dimensional vector £,■ in analogy to 
The components of e,- are or 1, and correspond to periodic or anti-periodic b.c, respectively. 
The partition function z e ' e2! for fixed £1,62,. . . can now be written as a linear combination of the 
partition functions Zjr f , e.g. for N = 1, 

Z ftA =4^-2j(-l)«S^-2^(-l)*5^+^(-l)« l ^^ iiS . (4.1) 

h h h,h 

As shown in [7], choosing periodic b.c. for all fermions in all directions, except antiperiodic in 
one direction for one single fermion, e.g. E\ = (0, 1) and £ (> i = (0,0), the corresponding partition 
function Z 0100 - vanishes at the massless, critical point, i.e. when the bare mass m is equal to the 
critical mass m c . Hence, the criterion can be used to determine m c for various couplings by tuning 
the bare mass m to the point where Z 0100 - = 0. 

In figure 2 we show the results of such a determination for the Thirring model. The plot on 
the left shows the partition function ratio Z 0100 /Z as a function of the bare mass am for square 
lattices with L = 8, . . . ,256 at the coupling g = 0.9. As the lattice size increases the jump from 
Z° 100 /Z ~ 1 to Z 0100 /Z ~ -0.25 is more and more pronounced, so that the critical point Z 0100 = 
can be determined very precisely. Since the critical point corresponds to a second order phase 
transition where the correlation length diverges, one should find a corresponding universal finite 
size scaling (FSS) behaviour. The inset in the left plot of figure 2 illustrates that this is indeed 
the case. There we show the partition function ratios as a function of (m — m c )L v ', and from the 
FSS we can determine the critical exponent v and the critical mass am c in the thermodynamic 
limit L — > 00 to a very high precision with a rather modest computational effort. The right plot in 



5 



Worm algorithm for the GN model 



Vidushi Maillart 




□ 


fermion 





boson 



<D a 
□ B 



Figure 3: In the left plot the single fermion and bound state correlation functions are plotted. The right plot 
shows the fermion and boson masses as a function of the bare mass am for the Thirring model at g = 1.1. 



figure 2 summarises our results for the critical mass in the thermodynamic limit as a function of 
the coupling g. 



5. Bound states of fermions 

It can easily be worked out that the correlation function of the pseudoscalar bound state Y75¥ 
in the Thirring model is given by 

(0{x)0{y)) = <1» 75 (y) 75 &60) • (5-1) 

In the loop formulation, this corresponds to two open fermion strings, one for each Majorana 
fermion flavour, with common endpoints at positions x and y. In complete analogy to the single 
fermion update using the fermionic 2-point correlation function, we can update the configurations 
using the bosonic bound state correlation function in eq.(5.1). In practice, we insert two instances 
of the bound state wave function &{x) = % j (x) y=, §2 (x) into the system and let them move around 
by employing again a local Metropolis update. This procedure samples the bound state correlation 
function, and at the same time updates the loop configuration of both fermion flavours. In order 
for this to work efficiently, the break-up and reconnection step described in section 3 is the crucial 
ingredient, since otherwise the algorithm would be restricted to move the bound state wave function 
only to sites where no fermion loop is present. Obviously, this would become increasingly difficult 
towards the critical point, where the fermion loops proliferate. 

The efficiency of the algorithm is illustrated in figure 3. In the left plot we show the single 
fermion and the bosonic bound state correlation functions at zero momentum obtained from a 
simulation of the Thirring model on a L = 128 lattice at coupling g = 1.1. It is remarkable that 
in both cases the signal can be followed over several orders of magnitude. Consequently, the 
corresponding masses can be reliably determined towards the continuum limit. This is illustrated 
in the right plot of figure 3 where we show the fermion and boson masses versus the bare mass am. 
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6. Conclusions and outlook 

Using Wilson's fermion discretisation, the path integral for the 0(2N) Gross-Neveu (GN) 
model can be described on the lattice in terms of interacting fermion loops. We discussed how the 
loop system can efficiently be simulated using open fermion strings. Single fermion and bound state 
correlation functions are measured while updating the system. In addition, the algorithm allows the 
direct calculation of ratios of partition functions with arbitrary fermion boundary conditions. We 
have successfully implemented the fermion loop algorithm for N = 1 , in which case the GN model 
is equivalent to the Thirring model, and presented first preliminary results for the determination of 
the critical point from the partition function ratios. Moreover, we also presented first promising 
results for the single fermion and the bound state masses. Currently we are working on measuring 
these quantities for the massive Thirring model in the continuum limit at various values of the 
couplings, in order to compare the results to predictions based on the equivalence of the model to 
the Sine-Gordon model. The extension of the algorithm to a larger number of fermions is interesting 
and rather straightforward. 

An obvious question to ask is whether and how the idea of updating an open fermionic string 
can be put to use in the context of gauged fermions or in higher dimensions. Successfull attempts 
were so far reported only in the strong coupling limit [1,2, 12], but there are many other interesting 
and promising extensions [13] using worm-type algorithms, even in connection with pure gauge 
theories [14]. 
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